#include "ClassMC.hpp"

double Likelihood_LSST_SNe( imcmc_double&   full_params,
                            double&         lndet,
                            double&         chisq,
                            void*           model,
                            void*           data,
                            istate&         state ){
    lndet = chisq = 0.0;

    CosmoTheory *cosmo = static_cast<CosmoTheory*>(model);
    LSST_SNe    *lsst = static_cast<LSST_SNe*>(data);

    double dmu;
    double *muzs = new double[lsst->A_ncol];
    double *muzp = new double[lsst->A_nrow];

    gsl_spline *spline_muzp = gsl_spline_alloc(gsl_interp_cspline, lsst->A_nrow);
    gsl_interp_accel *acc_muzp = gsl_interp_accel_alloc();

    for( int i=0; i<lsst->A_ncol; ++i ){
        muzs[i] = 5.0*log10( cosmo->engine->get_Dl(lsst->zs[i]) ) + 25.0;
    }

    for( int i=0; i<lsst->A_nrow; ++i ){
        muzp[i] = 0.0;
        for( int j=0; j<lsst->A_ncol; ++j ){
            muzp[i] += lsst->A(i,j)*muzs[j];
        }
    }

    gsl_spline_init(spline_muzp, lsst->zp, muzp, lsst->A_nrow);

    double MB = full_params["MB"];
    double MB_a = full_params["MB_a"];
    double MB_b = full_params["MB_b"];

    double mu, muerr;
    double MB_z;

    for( int i=0; i<lsst->sn_num; ++i ){

        MB_z = MB + lsst->sn_z[i]*(MB_a + MB_b*lsst->sn_z[i]) ;
        mu = gsl_spline_eval(spline_muzp, lsst->sn_z[i], acc_muzp);

        if( lsst->sn_z[i] <= lsst->muzperr_zmin )
            muerr = gsl_spline_eval(lsst->spline_muzperr, lsst->muzperr_zmin, lsst->acc_muzperr);
        else if( lsst->sn_z[i] >= lsst->muzperr_zmax )
            muerr = gsl_spline_eval(lsst->spline_muzperr, lsst->muzperr_zmax, lsst->acc_muzperr);
        else
            muerr = gsl_spline_eval(lsst->spline_muzperr, lsst->sn_z[i], lsst->acc_muzperr);

//      NOTE: the mock SNe data provides {zi,mu(zi)}, while the true observable is mb, the apparent magnitude.
//      Cosmological theory predicts only distance modulus, mu(z). We add to it the UNKNOWN absolute magnitude
//      MB(z), so we can predict mb(z).

        double mb = mu + MB_z;
        dmu = mb - ( lsst->sn_mu[i] + (-19.3) );
        chisq += pow(dmu/muerr,2);
    }

    delete[] muzs;
    delete[] muzp;
    gsl_spline_free(spline_muzp);
    gsl_interp_accel_free(acc_muzp);

    return -0.5*chisq;
}

/*
    JLA data has well measured redshifts, so that one does not need to consider
    making corrections to the predicted distance modulus, like what's needed for
    SNeIa samples with only photo-z.
*/
// double likelihood_SNeIa_JLA(    imcmc_double&   full_params,
//                                 double&         lndet,
//                                 double&         chisq,
//                                 void*           model,
//                                 void*           data,
//                                 istate&         state ){
//     lndet = chisq = 0.0;

//     CosmoTheory*    theory = static_cast<CosmoTheory*>(model);
//     SNeIa_JLA*      jla = static_cast<SNeIa_JLA*>(data);

//     double dmu;

//     for( int i=0; i<jla->sn_num; ++i ){
//         dmu = 5.0*log10( theory->engine->get_Dl(jla->sn_z[i]) ) + 25.0 - jla->sn_mu[i];
//         chisq += (dmu/jla->sn_err[i])*(dmu/jla->sn_err[i]);
//     }

//     return -0.5*chisq;
// }

// double likelihood_SNeIa(    imcmc_double&   full_params,
//                             double&         lndet,
//                             double&         chisq,
//                             void*           model,
//                             void*           data,
//                             istate&         state ){
//     lndet = chisq = 0.0;

//     CosmoTheory*    theory = static_cast<CosmoTheory*>(model);
//     SNeIa*          sndata = static_cast<SNeIa*>(data);

//     double chisq_temp, lndet_temp;

//     if( sndata->use_JLA ){
//         likelihood_SNeIa_JLA(full_params,lndet_temp,chisq_temp,theory,sndata->JLA,state);
//         lndet += lndet_temp;
//         chisq += chisq_temp;
//     }

//     if( sndata->use_LSST ){
//         likelihood_SNeIa_LSST(full_params,lndet_temp,chisq_temp,theory,sndata->LSST,state);
//         lndet += lndet_temp;
//         chisq += chisq_temp;
//     }

//     return -0.5*chisq;
// }
